1 Preface

The goal of this problem set is to develop some intuition about the impact of the number of nodes in the hidden layer of the neural network. We will use few simulated examples to have clear understanding of the structure of the data we are modeling and will assess how performance of the neural network model is impacted by the structure in the data and the setup of the network.

First of all, to compensate for lack of coverage on this topic in ISLR, let’s go over a couple of simple examples. We start with simulating a simple two class dataset in 2D predictor space with an outcome representative of an interaction between attributes. (Please notice that for the problems you will be working on this week you will be asked below to simulate a dataset using a different model.)

# fix seed so that narrative always matches the plots:
set.seed(1234567890)
nObs <- 1000
ctrPos <- 2
xyTmp <- matrix(rnorm(4*nObs), ncol = 2)
xyCtrsTmp <- matrix(sample(c( -1, 1)*ctrPos, nObs*4, replace = TRUE), ncol = 2)
xyTmp <- xyTmp + xyCtrsTmp
gTmp <- paste0("class", (1 + sign(apply(xyCtrsTmp, 1, prod)))/2)
plot(xyTmp, col = as.numeric(factor(gTmp)), pch = as.numeric(factor(gTmp)), xlab = "X1", ylab = "X2")
abline(h = 0)
abline(v = 0)

Symbol color and shape indicate class. Typical problem that will present a problem for any approach estimating a single linear decision boundary. We used similar simulated data for the random forest (week 10) problem set.

1.1 One hidden node

We can fit simple neural network (using all default values in the call to neuralnet – notice that both covariates and outcome have to be numeric as opposed to factor) and plot its layout (allowing for its output to be included in Rmarkdown generated report actually seems to be quite painful - one has to overwrite original implementation of plot.nn with the one that doesn’t call dev.new() that is included in this Rmarkdown file with echo = FALSE – to do the same you have to include that block into your Rmarkdown file also):

### Doesn't run: "requires numeric/complex ... arguments"
### nnRes <- neuralnet(g ~ X1 + X2, data.frame(g = gTmp, xyTmp))
nnRes <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)), xyTmp))
plot(nnRes)

That shows us a model with one node in a single hidden layer (default parameters).

We can lookup actual model predictions and recalculate them from input variables (in the field covariate) and model weight and activation function (fields weights and act.fct respectively):

head(nnRes$net.result[[1]])
##          [,1]
## [1,] 1.632381
## [2,] 1.632381
## [3,] 1.632377
## [4,] 1.547455
## [5,] 1.632381
## [6,] 1.632381
cbind(rep(1, 6), nnRes$act.fct(cbind(rep(1, 6), nnRes$covariate[1:6, ]) %*% nnRes$weights[[1]][[1]])) %*% nnRes$weights[[1]][[2]]
##          [,1]
## [1,] 1.632381
## [2,] 1.632381
## [3,] 1.632377
## [4,] 1.547455
## [5,] 1.632381
## [6,] 1.632381

Notice that input parameter linear.output governs whether activation function is called on the value of the output node or not:

nnResNLO <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)), xyTmp), linear.output = FALSE)
head(nnResNLO$net.result[[1]])
##           [,1]
## [1,] 0.9999904
## [2,] 0.9999904
## [3,] 0.9999904
## [4,] 0.9999904
## [5,] 0.9999903
## [6,] 0.9999904
nnResNLO$act.fct(cbind(rep(1, 6), nnResNLO$act.fct(cbind(rep(1, 6), nnResNLO$covariate[1:6, ]) %*% nnResNLO$weights[[1]][[1]])) %*% nnResNLO$weights[[1]][[2]])
##           [,1]
## [1,] 0.9999904
## [2,] 0.9999904
## [3,] 0.9999904
## [4,] 0.9999904
## [5,] 0.9999903
## [6,] 0.9999904
quantile(nnResNLO$net.result[[1]])
##        0%       25%       50%       75%      100% 
## 0.9999903 0.9999904 0.9999904 0.9999904 0.9999904

As the last statement (the quantiles of the predicted out) above shows, the use of activation function limiting predicted values to \([0;1]\) range when modeling outcome taking values outside of \([0;1]\) interval does not result in a very useful model. In this case with true outcome values constrained to \(\{1, 2\}\) so that the error is minimized by predicting every outcome to be as close to \(1\) as possible.

Using binary – 0 or 1 – outcome produces more useful model when activation function is applied to the output node (linear.output = FALSE) and allows use of cross-entropy error function (often used in classification setting in combination with the activation function applied to the output layer):

nnResNLO01 <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)) - 1, xyTmp), linear.output = FALSE)
quantile(nnResNLO01$net.result[[1]], c(0, 0.1, 0.25, 0.5, 1))
##         0%        10%        25%        50%       100% 
## 0.02967696 0.03032911 0.57518638 0.63188940 0.63190537
if ( FALSE ) {
  # from v.3.5.x(?) this started to throw an error instead of silent switch to "sse":
  nnResNLO12CE <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)), xyTmp), linear.output = FALSE, err.fct = "ce")
}
nnResNLO01CE <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)) - 1, xyTmp), linear.output = FALSE, err.fct = "ce")
head(nnResNLO01CE$net.result[[1]])
##           [,1]
## [1,] 0.3722697
## [2,] 0.3773444
## [3,] 0.9852998
## [4,] 0.6211743
## [5,] 0.3722790
## [6,] 0.3722697
nnResNLO01CE$act.fct(cbind(rep(1, 6), nnResNLO01CE$act.fct(cbind(rep(1, 6), nnResNLO01CE$covariate[1:6, ]) %*% nnResNLO01CE$weights[[1]][[1]])) %*% nnResNLO01CE$weights[[1]][[2]])
##           [,1]
## [1,] 0.3722697
## [2,] 0.3773444
## [3,] 0.9852998
## [4,] 0.6211743
## [5,] 0.3722790
## [6,] 0.3722697
quantile(nnResNLO01CE$net.result[[1]])
##        0%       25%       50%       75%      100% 
## 0.3722697 0.3722697 0.3722860 0.4137822 0.9853637

We can plot model output indicating class identity (left panel below) that tells us that when true outcome values are constrained to 1 or 2, sum of squared errors is used as error function and output node values are used as-is (not transformed by activation function – linear.output = TRUE by default), the majority of the points were estimated to be close to 1 or 1.6 and that majority of those estimated to be close to 1 correspond to the about half of the observations at the first level of the class factor (i.e. numerical value of 1). It is also easy to see that those with predicted value of 1.6 represent roughly 1:2 mix of observations from the first and second levels of the outcome respectively, so that \(1.6 \approx (1 + 2*2)/3 = 5/3\) approximately equals average of their numerical values corresponding to the levels of the factor representing them.

The nature of the model estimated by neuralnet in this (very simple!) case becomes even more intuitive if we render all points in the area encompassing our training set with model predictions and overlay training dataset on top of that (right panel below). It is immediately apparent that this model identified a line in this 2D space separating one cloud of points belonging mostly to one class from all others so that predicted values are approximately equal to the average outcome on each side of this decision boundary.

plotNNpreds2D2class <- function(inpNN, inpClassThresh, inpGrid = ( -60:60)/10) {
  tmpClrPch <- as.numeric(factor(inpNN$response))
  plot(inpNN$net.result[[1]], col = tmpClrPch, pch = tmpClrPch)
  table(inpNN$net.result[[1]][, 1]>inpClassThresh, inpNN$response)
  xyGridTmp <- cbind(X1 = rep(inpGrid, length(inpGrid)), X2 = sort(rep(inpGrid, length(inpGrid))))
  # before predict.nn existed:
  #gridValsTmp <- compute(inpNN, xyGridTmp)$net.result
  gridValsTmp <- predict(inpNN, xyGridTmp)
  errTmp <- sum(inpNN$err.fct(inpNN$net.result[[1]][, 1], inpNN$response))
  plot(xyGridTmp, col = as.numeric(gridValsTmp>inpClassThresh) + 1, pch = 20, cex = 0.3, main = paste("Error:", round(errTmp, 6)))
  points(inpNN$covariate, col = tmpClrPch, pch = tmpClrPch)
  ## Equations defining decision boundary:
  ## 1*w0 + X1*w1 + X2*w2 = 0, i.e.:
  ## 0 = inpNN$weights[[1]][1] + inpNN$weights[[1]][2]*X1 + inpNN$weights[[1]][3]*X2, i.e:
  ## X2 = ( - inpNN$weights[[1]][1] - inpNN$weights[[1]][2]*X1) / inpNN$weights[[1]][3]
  for ( iTmp in 1:ncol(inpNN$weights[[1]][[1]]) ) {
    abline( - inpNN$weights[[1]][[1]][1, iTmp] / inpNN$weights[[1]][[1]][3, iTmp], -inpNN$weights[[1]][[1]][2, iTmp] /inpNN$weights[[1]][[1]][3, iTmp], lty = 2, lwd = 2)
  }
}
old.par <- par(mfrow = c(1, 2), ps = 16)
plotNNpreds2D2class(nnRes, 1.3)

par(old.par)

Similar, aside from a different position of decision boundary, results can be obtained when using binary representation of the outcome with cross-entropy as error function together with applying activation function to the output layer:

plot(nnResNLO01CE)

Once activation function is applied to the output node, the output values are bound to the \([0, 1]\) interval. The predicted values that are far enough from the decision boundary are also approximately set to the average of the outcome in that subspace:

old.par <- par(mfrow = c(1, 2), ps = 16)
plotNNpreds2D2class(nnResNLO01CE, 0.5)

par(old.par)

The important points resulting from the results shown in the figures above are the following:

  • as simple of a model as the one that was employed here (with one node in a single hidden layer along with all other default parameters) cannot do much better than what we observed here
  • because calling (default-logistic) activation function on a given linear combination of the input variables more or less amounts to assigning almost all points on one side of hyperplane (line in 2D, plane in 3D, etc.) to zero and on the other side – to unity
  • weights involved in transforming outcome of the hidden layer into model predictions will change those zeroes and ones to values closer to the desired outcome values, but still, use of such a simple model (with a single hidden node) employed here to prime our intuition more or less amounts to splitting covariate space into two half spaces by a hyperplane and assigning almost constant outcomes to the vast majority of the points on either side of it
  • the weights for the inputs into the single hidden node shown in the network layout plot above and stored in weights field of the result returned by neuralnet define this hyperplane (line in 2D, etc.) shown as dashes in the panels on the right above
  • this hyperplane is where sum of weighted input variables and an intercept is identical zero (and thus the result of logistic activation function is 0.5 rapidly becoming zero or one for points further away from this boundary)

1.2 Two hidden nodes, single hidden layer

Now, let’s add another node to the hidden layer of this network. From the above, we know what to expect as a result of that – another hyperplane (line in 2D) will be added to the space of covariates now dividing it into (depending on whether those hyperplans are almost parallel or not) three or four subspaces, consequently, assigning most of the points to three or four potentially different constants. Clearly, this level of granularity could suffice for developing a model that would do quite well in our toy example.

To have more than one node in the hidden layer we set hidden parameter to the number of nodes in it (length of vector provided as hidden parameter governs the number of the hidden layers in the network – we still use one layer here):

set.seed(1234567)
nnRes2 <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)), xyTmp), hidden = 2)
plot(nnRes2)

We can see that now resulting network has two nodes in a single hidden layer, which two covariates enter with weights that are approximately comparable in magnitude and opposite in sign. Their comparably weighted sum added to a constant close to one gives the outcome value of this model. The effect of those weights in defining decision boundaries in the space of predictors is best seen from the figure below:

old.par <- par(mfrow = c(1, 2), ps = 16)
plotNNpreds2D2class(nnRes2, 1.5)

par(old.par)

This model sets up two almost parallel lines that encompass most of the observations from the first class, leaving most of the observations from the second class outside of the resulting slab. Now let’s repeat fitting neural network three times (each time starting with random choice of starting weights in the model) and compare stability of the resulting models:

old.par <- par(mfcol = c(2, 3), ps = 16)
for ( iTry in 1:3 ) {
  nnRes2 <- neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)), xyTmp), hidden = 2)
  plotNNpreds2D2class(nnRes2, 1.5)
}

par(old.par)

We can see that quite frequently given the parameters used the process converges to a suboptimal solution with about half of the observations remaining in “gray” zone where their assignment to either of the classes is not immediately apparent.

Aside from the multitude of local minima for neural network fitting procedure that could prevent it from finding better solutions, the main point to take from this exercise is that adding more nodes to the hidden layer (with all other default choices employed here) amounts to adding more hyperplanes bisecting the space of predictors, creating more and more subspaces where the outcome can take different values (often close to a constant in each subspace). Obviously, the geometry of the resulting decision surfaces can become quite complicated even with modest number of nodes in the hidden layer. Lastly, these considerations provide some intuition for considering what could be a useful number of hidden nodes in the model. In thinking about that it might be useful to consider how many such hyperplanes could be sufficient to effectively separate observations belonging to different outcome categories. Not that we necessarily would have such knowledge ahead of time, but this might prove to be a complementary way to think about the problem in addition to the often sited empirical guidelines that are based on the number of predictor variables, etc.

The point of this problem set is to assess how these aspects of neural network fitting play out in another simulated dataset.

Lastly, we also would like you to consider what combination of error function and output node transformation you would like to use for this week problem set. Below are three calls of neuralnet timed using different combinations of err.fct and linear.output:

system.time(invisible(neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)) - 1, xyTmp), hidden = 2, linear.output = TRUE, err.fct = "sse")))
##    user  system elapsed 
##   0.227   0.017   0.244
system.time(invisible(neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)) - 1, xyTmp), hidden = 2, linear.output = FALSE, err.fct = "sse")))
##    user  system elapsed 
##   1.560   0.174   1.734
system.time(invisible(neuralnet(g ~ X1 + X2, data.frame(g = as.numeric(factor(gTmp)) - 1, xyTmp), hidden = 2, linear.output = FALSE, err.fct = "ce")))
## Warning: Algorithm did not converge in 1 of 1 repetition(s) within the
## stepmax.
##    user  system elapsed 
##  17.526   1.893  19.424

Obviously, for our toy example use of untransformed outcome from the output node and sum of squares as error function results in the fastest convergence. But then the decision as to how to translate its predictions to class categories is yours. On the other hand, one might argue that more principled approach since we are dealing with classification problem would be to use logistic transform of the outcome to contain it within \([0;1]\) interval and cross-entropy as an error function. Except that it seems to run nearly two orders of magnitude slower in this case and does not necessarily converge with default values for the rest of the arguments. Please feel free to experiment with how they will perform for the problems you are presented with below and decide for yourself which combination of these parameters is more suitable for the task at hand.

2 Problem 1 (10 points): generate 3D data with nested cubes

Simulate data with n = 1000 observations and p = 3 covariates – all random variables from uniform distribution contained on \([ -1, 1]\) interval: \(X_i\sim \mathcal{U}( -1, 1), i = \{1, 2, 3\}\). Resulting observations should be uniformly filling in a 3D cube centered at \((0, 0, 0)\) with the edge of length \(2\). Create two category class variable assigning all observations within a smaller cube – also centered at \((0, 0, 0)\) with the edge length of \(2^\frac{2}{3}\) – to one class category and observations outside of this cube – to the second class category.

Please note that this dataset is entirely different from the one used in the preface – you will need to write code simulating it on your own. Somewhat related 2D dataset was used as a motivational example at week 12 (SVM) lecture before introducing kernels and SVMs. However, the example in the lecture was in 2D (three-dimensional problem here), had spheric boundary (here we work with nested cubes) and used normal (here uniform) distribution. Since you will be reusing this code in the following two problems it is probably best to turn this procedure into a function with appropriate parameters.

Ten points available for this problem are composed of accomplishing the following tasks:

  1. Correct implementation of the function generating data as prescribed above (2 points)
library(scatterplot3d)
set.seed (1)
nObs <- 1000
x <- runif(nObs, -1, 1)
y <- runif(nObs, -1, 1)
z <- runif(nObs, -1, 1)

# The inner cube is 
myBound <- 2^(2/3) * 0.5

twoClassesIn3D <- cbind(x,y,z, (x > myBound | y > myBound | z > myBound | x < -myBound | y < -myBound | z < -myBound ));
  1. Check and demonstrate numerically that the resulting class assignment splits these observations, subject to sampling variability, evenly between these two groups (2 points)
  2. Explain why this is achieved by using \(2^\frac{2}{3}\) as the length of the inner cube edge that defines true decision boundary (2 points)
print("Count of points in class 1 divided by count of points in the other class:")
## [1] "Count of points in class 1 divided by count of points in the other class:"
sum(twoClassesIn3D[,4]) / length(twoClassesIn3D[,4])
## [1] 0.513
print("Volume of inner cube = the cube of the shape's edge length:")
## [1] "Volume of inner cube = the cube of the shape's edge length:"
(2^(2/3))^3
## [1] 4
print("Volume of the remaining space is the total outer cube volume (edge=2) minus the inner cube volume (edge=2^(2/3)):")
## [1] "Volume of the remaining space is the total outer cube volume (edge=2) minus the inner cube volume (edge=2^(2/3)):"
(2^3) - (2^(2/3))^3
## [1] 4

As we can see in the above result, classes are split about 50-50. This is because the volume of the two spaces is equal (4 cubic units).

  1. Plot values of the resulting covariates projected at each pair of the axes indicating classes to which observations belong with symbol color and/or shape (you can use function pairs, for example) (2 points)
scatterplot3d(twoClassesIn3D[,1],
              twoClassesIn3D[,2],
              twoClassesIn3D[,3],
              xlab = "x",
              ylab = "y",
              zlab = "z",
              main = "Two classes uniformly distributed in three dimensions",
              
              # FIXME: LaTeX support in R plots?
              #sub = "red points are inside the $2^\frac{2}{3}$ cube at origin",
              
              sub = "Red points are inside the cube of length 2^0.6666 at origin",
              color = ifelse(twoClassesIn3D[,4], 'red', 'black'))

  1. Reflect on the geometry of this problem by answering the following question: what is the smallest number of planes in 3D space that would completely enclose points from the “inner” class?

Four planes – four triangular faces create an enclosure.

Is this number equal to the number of cube faces or is it something smaller? Larger? Please note that “enclose” above does not mean “perfectly discriminate between the points assigned to two classes” (2 points)

Smaller than the face count of a cube. Cube has six faces.

3 Problem 2 (20 points): neural network classifier

For the dataset simulated above fit neural networks with 1 through 6 nodes in a single hidden layer (use neuralnet implementation). For each of them calculate training error (see an example in Preface where it was calculated using err.fct field in the result returned by neuralnet).

# Get training and test data from wifiLocDat list
getTestData <- function(nObsArg = 1000) {

  nObs <- nObsArg
  x <- runif(nObs, -1, 1)
  y <- runif(nObs, -1, 1)
  z <- runif(nObs, -1, 1)
  
  # The inner cube boundary
  myBound <- 2^(2/3) * 0.5
 
  twoClassesIn3D <- cbind(x,y,z, (x > myBound | y > myBound | z > myBound | x < -myBound | y < -myBound | z < -myBound ));

  colnames(twoClassesIn3D) <- c("x", "y", "z", "outcome")
  
  # Get random indices
  randomIndices <- sample(1:nObs, floor(nObs/2))

  # Get indices not in randomIndices
  randomIndicesInverse <- setdiff(seq(1:nObs), randomIndices)

  xyzTrain   <- twoClassesIn3D[ randomIndices, 1:3]
  xyzTest    <- twoClassesIn3D[-randomIndices, 1:3]
  classTrain <- twoClassesIn3D[ randomIndices, 4]
  classTest  <- twoClassesIn3D[-randomIndices, 4]

  return ( list(
    "xyzTest"    = xyzTest,
    "xyzTrain"   = xyzTrain,
    "classTest"  = classTest,
    "classTrain" = classTrain
  ))
}

gtd <- getTestData()
xyzTest    = gtd$xyzTest
xyzTrain   = gtd$xyzTrain
classTest  = gtd$classTest
classTrain = gtd$classTrain
computeErrorRate <- function (myTable) {
  return (1 - sum(diag(myTable)) / sum(myTable))
}

getNnTmpTable <- function(xyzTest, xyzTrain, classTest, classTrain, hiddenArg) {
  nnRes <- neuralnet(outcome ~ ., cbind(xyzTrain, outcome = classTrain), hidden = c(hiddenArg), rep = 1, stepmax=1e+06, threshold = 0.1, linear.output = FALSE, err.fct = 'ce')
  nnTmpTbl <- table(classTest, round(predict(nnRes, newdata = xyzTest)))
  return (nnTmpTbl)
}

myHeights <- c()
myHiddenLayerNodeCounts <- 1:6

for (hA in myHiddenLayerNodeCounts) {
  nnTable <- getNnTmpTable(xyzTest, xyzTrain, classTest, classTrain, hiddenArg = hA)
  print(paste("Number of nodes in hidden layer: ", hA))
  
  newHeight <- computeErrorRate(nnTable)
  print(paste("Error rate: ", newHeight))
  myHeights <- c(myHeights, newHeight)
}
## [1] "Number of nodes in hidden layer:  1"
## [1] "Error rate:  0.416"
## [1] "Number of nodes in hidden layer:  2"
## [1] "Error rate:  0.3"
## [1] "Number of nodes in hidden layer:  3"
## [1] "Error rate:  0.24"
## [1] "Number of nodes in hidden layer:  4"
## [1] "Error rate:  0.084"
## [1] "Number of nodes in hidden layer:  5"
## [1] "Error rate:  0.204"
## Warning: Algorithm did not converge in 1 of 1 repetition(s) within the
## stepmax.
## Error in cbind(1, pred) %*% weights[[num_hidden_layers + 1]]: requires numeric/complex matrix/vector arguments
# Plot the error rates for each hidden layer node count
barplot(height = myHeights,
        names.arg = myHiddenLayerNodeCounts,
        col = rainbow(length(myHeights)),
        xlab = "Number of nodes in hidden layer",
        ylab = "Error rate")
## Error in barplot.default(height = myHeights, names.arg = myHiddenLayerNodeCounts, : incorrect number of names

# Classification (so linear.output = TRUE), not regression
# for (hiddenLayerNode in 1:6) {
# for (hiddenLayerNode in 1:6) {
#   nnRes <- neuralnet(outcome ~ ., twoClassesIn3D, hidden = c(hiddenLayerNode), linear.output = FALSE)
#   print(nnRes)
#   plot(nnRes)
# }

Simulate another independent dataset (with n = 10,000 observations to make resulting test error estimates less variable) using the same procedure as above (3D, two classes as nested cubes) and use it to calculate the test error at each number of hidden nodes. Plot training and test errors as function of the number of nodes in the hidden layer. What does resulting plot tells you about the interplay between model error, model complexity and problem geometry? What is the geometrical interpretation of this error behavior?

gtd <- getTestData(10000)
xyzTest    = gtd$xyzTest
xyzTrain   = gtd$xyzTrain
classTest  = gtd$classTest
classTrain = gtd$classTrain

for (hA in 1:6) {
  nnTable <- getNnTmpTable(xyzTest, xyzTrain, classTest, classTrain, hiddenArg = hA)
  print(paste("Number of nodes in hidden layer: ", hA))
  print(paste("Error rate: ", computeErrorRate(nnTable)))
}
## [1] "Number of nodes in hidden layer:  1"
## [1] "Error rate:  0.4034"
## [1] "Number of nodes in hidden layer:  2"
## [1] "Error rate:  0.3008"
## [1] "Number of nodes in hidden layer:  3"
## [1] "Error rate:  0.3032"
## [1] "Number of nodes in hidden layer:  4"
## [1] "Error rate:  0.1562"
## Warning: Algorithm did not converge in 1 of 1 repetition(s) within the
## stepmax.
## Error in cbind(1, pred) %*% weights[[num_hidden_layers + 1]]: requires numeric/complex matrix/vector arguments

4 Problem 3 (30 points): evaluate impacts of sample size and noise

Setup a simulation repeating procedure described above for n = 100, 200 and 500 observations in the training set as well as adding none, 1, 2 and 5 null variables to the training and test data (and to the covariates in formula provided to neuralnet). Draw values for null variables from uniform distribution on \([ -1, 1]\) as well and do not use them in the assignment of the observations to the class category (e.g. x <- matrix(runif(600, -1, 1), ncol = 6); cl <- as.numeric(factor(rowSums(abs(x[, 1:3])<2^( -1/3)) = =3)) creates dataset with three informative and three null variables). Repeat calculation of training and test errors at least several times for each combination of sample size, number of null variables and size of the hidden layer simulating new training and test dataset every time to assess variability in those estimates. Present resulting error rates so that the effects of sample size and fraction of null variables can be discerned and discuss their impact of the resulting model fits.

5 Extra 10 points problem: model WiFi localization data

Use neuralnet to model the outcome in WiFi localization dataset that we used in previous weeks using outcome in its original, four-levels format. Your neuralnet models will need to have four output nodes to represent outcome with four levels. Obtain training and test error for this model for several sizes of the hidden layer. Compare resulting test error in predicting Location 3 to that observed for random forest, SVM and KNN approaches earlier in the course.

Session info

For reproducibility purposes it is always a good idea to capture the state of the environment that was used to generate the results:

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.14.2
## 
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] scatterplot3d_0.3-41 ggplot2_3.2.1        neuralnet_1.44.2    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2       knitr_1.24       magrittr_1.5     tidyselect_0.2.5
##  [5] munsell_0.5.0    colorspace_1.4-1 R6_2.4.0         rlang_0.4.1     
##  [9] stringr_1.4.0    dplyr_0.8.3      tools_3.5.1      gtable_0.3.0    
## [13] xfun_0.9         withr_2.1.2      htmltools_0.4.0  assertthat_0.2.1
## [17] yaml_2.2.0       lazyeval_0.2.2   digest_0.6.22    tibble_2.1.3    
## [21] crayon_1.3.4     purrr_0.3.2      glue_1.3.1       evaluate_0.14   
## [25] rmarkdown_1.15   stringi_1.4.3    compiler_3.5.1   pillar_1.4.2    
## [29] scales_1.0.0     pkgconfig_2.0.3

The time it took to knit this file from beginning to end is about (seconds):

proc.time() - ptStart
##     user   system  elapsed 
## 1697.848  732.501 2196.329